## ----------------------------------------------------------------------
##
## Regress ninka-dantai outcomes on raid damage
## Note:
##
## ----------------------------------------------------------------------
## environment setting
Sys.setenv(LANGUAGE="en")
gc();gc()
rm(list = ls())
options(scipen = 999)   ## Disable scientific notation

## ----------------------------------------------------------------------
##
## Initial settings
##
## ----------------------------------------------------------------------
## Set directory
## ----------------------------------------------------------------------
root_dir = "YOUR_DIRECTORY/Replication_Folder"
## Load packages and functions
source(file.path(root_dir, "2_Code/1_PackagesFunctions.R"))

## ----------------------------------------------------------------------
## Data directory
data_dir = root_dir	%>%	file.path("1_Data")	%>%	dir_create()
## Output directory
output_dir = root_dir	%>%	file.path("3_Result")	%>%	dir_create()
figure_dir = root_dir	%>%	file.path("4_Figures")	%>%	dir_create()
## Working directory
working_dir = root_dir	%>%	file.path("5_Working")	%>%	dir_create()
## ----------------------------------------------------------------------


## ----------------------------------------------------------------------
## Load data
## ----------------------------------------------------------------------
## --- Census data
tokyo_census_tbl = data_dir	%>%	file.path("census_data4.csv")	%>%
  ##census_data4
  read_csv()	%>%	
  filter(year == 2015) %>%
  mutate(key_code = as.character(key_code))

## --- All years
tokyo_ninka_tbl = data_dir	%>%	file.path("test_ninka_2.dta")	%>%	read_dta()	%>%
  rename(
    dummy_ninka = d_ninka,
    ln_ninka_days = lnninka_days
  )
## --- Since 1993
tokyo_ninka_tbl_1993 = data_dir	%>%	file.path("test_ninka_sinceFY1993_2.dta")	%>%	read_dta()	%>%
  rename(
    dummy_ninka = d_ninka,
    ln_ninka_days = lnninka_days
  )
## --- Sports
# sports_tbl = root_dir	%>%	file.path("@@SUBMISSION/JOP_RR1/supporting_materials/R2_Q5a/sports_club_list.csv")	%>%
sports_tbl = data_dir	%>%	file.path("sports_club_list.csv")	%>%
  read_csv(col_types = cols())	%>%
  rename(
    club_name = 1, cityName = ku_name,
    chochoazaName = address_jp)	%>%
  mutate(
    chochoazaName = if_else(
      str_detect(chochoazaName, pattern = "\\d"),
      true = str_c(chochoazaName, "丁目"),
      false = chochoazaName),
    chochoazaName = stringi::stri_trans_general(chochoazaName, id = "Halfwidth-Fullwidth"),
    SportClubDummy = 1,
    zipcode = str_remove_all(zipcode, pattern = "-")
  )

## ----------------------------------------------------------------------
##
## Prepare regression-related objects
##
## ----------------------------------------------------------------------
## Link function
binom_link = binomial(link = logit)
count_link = quasipoisson
## ----------------------------------------------------------------------
## Outcome and treatment variables
## ----------------------------------------------------------------------
depvar_vec_combined = c("ninka_days", "ln_ninka_days", "dummy_ninka", "n_ninka")
placebo_depvar_vec = c("SportClubDummy", "LocalSportClubDummy")
treatment_vec = c("ratio_damage", "ln_damage_ratio", "binary_damage_ratio")
robust_treatment_vec = c(treatment_vec[2:3], str_c("as.factor(", treatment_vec[1],")"))
## ----------------------------------------------------------------------
## Covariates
## ----------------------------------------------------------------------
## District FE
fe_term = "DistrictFE"
## Geographical features
spatial_covariates = c(
  ## Residential ratio and geographical area
  "ratio_residential", "poly_area_ln", "n_neighbor",
  ## Prewar population density
  "PopDensity_1939_km2_ln",
  ## Terrain variables
  "mean_elevation_ln", "mean_slope_ln", "river_distance_ln",
  ## Railway and train stations
  "railway_length_ln", "n_stations", "SpLag_railway_length_ln", "SpLag_n_stations"
)
## ----------------------------------------------------------------------
## Location variables
distance_termz_base = c("palace_dist_ln", "minTargetDistance_ln")
polynomial_degree = 5
## --- Polynomials and splines (single dimension)
distance_term_specification = list(
  ## --- Polynomials
  str_c("poly(", distance_termz_base, ", degree = ", polynomial_degree, ", simple = TRUE, raw = TRUE)"),
  ## --- Splines
  str_c("s(", distance_termz_base, ")"),
  ## --- none (comparison purposes)
  "",
  ## --- 3rd order polynomials
  str_c("poly(", distance_termz_base, ", degree = 3, simple = TRUE, raw = TRUE)")
)
## --- Lon-lat (two-dimensional) term
spatial_term_specification = c(
  ## --- Polynomials
  str_c("poly(z_lon, z_lat, degree = ",  polynomial_degree, ", simple = TRUE, raw = TRUE)"),
  ## --- Splines
  "te(longitude, latitude)",
  ## --- none (comparison purposes)
  "",
  ## --- 3rd order polynomials
  "poly(z_lon, z_lat, degree = 3, simple = TRUE, raw = TRUE)"
)
## ----------------------------------------------------------------------
## Extra (two-dimensional) polynomials
extra_covariates_specification = c(
  ## --- none (baseline)
  "",
  ## --- with extra terms (robustness check)
  str_c("poly(ratio_residential, PopDensity_1939_km2_ln, degree = ", polynomial_degree, ", simple = TRUE, raw = FALSE)"))
## ----------------------------------------------------------------------

## ----------------------------------------------------------------------
## Combine the post-1993 records and full records
## ----------------------------------------------------------------------
depvar_1993_tbl = tokyo_ninka_tbl_1993	%>%
  select(
    key_code, chochoazaName, all_of(depvar_vec_combined)
  )	%>%
  rename_at(vars(ninka_days:n_ninka), ~ str_c(., "_1993"))
## ----------------------------------------------------------------------
tmp_combined_tbl = tokyo_census_tbl	%>%
  left_join(
    tokyo_ninka_tbl	%>%
      select(key_code, chochoazaName, all_of(depvar_vec_combined))
  )	%>%
  left_join(depvar_1993_tbl, by = c("key_code", "chochoazaName"))	%>%
  left_join(sports_tbl, by = c("cityName", "chochoazaName", "zipcode"))	%>%
  # left_join(full_poly, by = "key_code")	%>%
  rownames_to_column(var = "row_count")	%>%
  mutate(
    DistrictFE = as.factor(prewar_district),
    SportClubDummy = if_else(is.na(SportClubDummy), true = 0, false = SportClubDummy)
  )	%>%
  select(row_count, everything())
## ----------------------------------------------------------------------

## ----------------------------------------------------------------------
## Spatially-lagged variables
## ----------------------------------------------------------------------
full_poly = data_dir	%>%	file.path("RaidShp_May2020b.rds")	%>%
  read_rds()	%>%
  filter(row_id!=88) %>% ## drop palace polygon
  select(key_code, ratio_damage)	%>%
  left_join(tmp_combined_tbl	%>%	select(key_code, SportClubDummy))	%>%
  mutate(
    SportClubDummy = if_else(is.na(SportClubDummy), true = 0, false = SportClubDummy)
  )
## Use different SWM
full_contig_nb = poly2nb(full_poly, queen = TRUE)
full_contig_listw = nb2listw(full_contig_nb, style = "W", zero.policy = TRUE)	## 1 neighborhood without neighbors
## Lagged damage SWM
full_contig_nb_damage = poly2nb(full_poly	%>%	drop_na(ratio_damage), queen = TRUE)
full_contig_listw_damage = nb2listw(full_contig_nb_damage, style = "W", zero.policy = TRUE)	## 1 neighborhood without neighbors
# lapply(full_contig_listw_damage$weights, length)	%>%	unlist()	## N neighbor cells
## ----------------------------------------------------------------------
## Lagged damage
lagged_damage_tbl = full_poly	%>%	drop_na(ratio_damage)	%>%
  mutate(
    LaggedDamage = lag.listw(full_contig_listw_damage, var = full_poly	%>%	drop_na(ratio_damage)	%>%	pull(ratio_damage),
                             zero.policy = TRUE)
  )	%>%
  as_tibble()	%>%
  select(key_code, LaggedDamage)
## Construct lagged outcome
lagged_variables_tbl = full_poly	%>%
  mutate(
    LocalSportClubDummy = lag.listw(full_contig_listw, var = full_poly$SportClubDummy, zero.policy = TRUE, NAOK = TRUE),
    LocalSportClubDummy = if_else(SportClubDummy + LocalSportClubDummy > 0, true = 1, false = 0)
  )	%>%
  as_tibble()	%>%
  select(key_code, SportClubDummy, LocalSportClubDummy)	%>%
  left_join(lagged_damage_tbl, by = "key_code")
lagged_variables_tbl	%>%	summary()
## ----------------------------------------------------------------------

## ----------------------------------------------------------------------
## Combine and save
## ----------------------------------------------------------------------
combined_tbl = tmp_combined_tbl	%>%	left_join(lagged_variables_tbl)	%>%
  mutate(
    zipcode = str_c(str_sub(zipcode, start = 1, end = 3), "-", str_sub(zipcode, start = 4, end = 7)),
    zipcode = if_else(zipcode == "962-0838", true = "173-0027", false = zipcode)
  )	%>%
  select(
    row_id, row_count:pop, zipcode, -year,
    all_of(depvar_vec_combined), contains("_1993"),
    all_of(placebo_depvar_vec), all_of(treatment_vec), LaggedDamage,
    all_of(fe_term), all_of(spatial_covariates), all_of(distance_termz_base),
    z_lon, z_lat, longitude, latitude,
    palace_dist)	%>%
  drop_na(ratio_damage)	%>%
  write_csv(file = data_dir	%>%	file.path("ANA_data.csv"))

## Read zipcode-level data
zipcode_tbl = data_dir	%>%	file.path("ANA_zipcode_data.csv")	%>%	read_csv(col_types = cols())

## ----------------------------------------------------------------------



## ----------------------------------------------------------------------
## Baseline formula objects
## ----------------------------------------------------------------------
mdl_lst = list(
  ## NOT IN USE, Polynomial: lm()-compatible specification
  polynomial_base = chr2fml(
    "outcome",
    idv_list = list(
      treatment_vec[1],
      spatial_covariates, fe_term,
      distance_term_specification[[1]],
      spatial_term_specification[1],
      extra_covariates_specification[1]	## none
    )
  ),
  ## Polynomial: felm()-compatible specification
  polynomial_base_felm = chr2fml_felm(
    "outcome",
    idv_list = list(
      treatment_vec[1],
      spatial_covariates,
      distance_term_specification[[1]],
      spatial_term_specification[1],
      extra_covariates_specification[1]	## none
    ),
    fe_list = fe_term,
    se_cluster = "row_count"	## Robust SE
    # se_cluster = fe_term	## Cluster SE
  ),
  ## Spline: gam()-compatible specification
  gam_base = chr2fml(
    "outcome",
    idv_list = list(
      treatment_vec[1],
      spatial_covariates, fe_term,
      distance_term_specification[[2]],
      spatial_term_specification[2],
      extra_covariates_specification[1]	## none
    )
  ),
  ## gam with spatial 5th-order poly (for tab a9)
  gam_sp5th = chr2fml(
    "outcome",
    idv_list = list(
      treatment_vec[1],
      spatial_covariates, fe_term,
      distance_term_specification[[1]],
      spatial_term_specification[1],
      extra_covariates_specification[1]	## none
    )
  )
)
## ----------------------------------------------------------------------


## ----------------------------------------------------------------------
##
## Table 1
## Result I: Baseline regressions
##
## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[3], " ~ ."))	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[3], "_1993 ~ ."))	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[2], " ~ ."))	%>%
    # update(str_c(depvar_vec_combined[1], " ~ ."))	%>%
    gam(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[4], " ~ ."))	%>%
    gam(data = combined_tbl, family = count_link)
)

## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table 1: Panel A
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)

## ----------------------------------------------------------------------
## Panel B: Dummy regressions
## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst_dummy = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(
        depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3]
      )	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(
        depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3]
      )
    )	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(
      str_c(
        depvar_vec_combined[2], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3]
      )
    )	%>%
    # update(str_c(depvar_vec_combined[1], " ~ ."))	%>%
    gam(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(
        depvar_vec_combined[4], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3]
      )
    )	%>%
    gam(data = combined_tbl, family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table 1: Panel B
stargazer(
  rslt_lst_dummy,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------




## ----------------------------------------------------------------------
##
## Table A2
## Descriptive Statistics (1 of 2)
## See 3_CensusMainTxt&AppxA for the descriptive stats for
## the outcome variables created from census.
##
## ----------------------------------------------------------------------
## Summary statistics
filtered_spatial_covariates <- spatial_covariates %>% 
  as.data.frame() %>% 
  filter(. != "ratio_residential") %>% 
  t() %>% 
  as.vector()
tableContinuous(
  vars = combined_tbl	%>%
    select(
      all_of(treatment_vec[1]), LaggedDamage,
      ln_ninka_days, dummy_ninka, dummy_ninka_1993, n_ninka,
      all_of(placebo_depvar_vec),
      ratio_residential, all_of(distance_termz_base),
      all_of(filtered_spatial_covariates), 
      longitude, latitude)	%>%	as.data.frame(),
  group = combined_tbl$ratio_damage > median(combined_tbl$ratio_damage),
  stats = c("n", "mean", "s", "median", "iqr"),
  prec = 3)
## ----------------------------------------------------------------------



## ----------------------------------------------------------------------
##
## Table A4
## Result V: >median damage dummy
##
## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ratio_damage + binary_damage_ratio")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[3], "~ . - ratio_damage + binary_damage_ratio"))	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ . - ratio_damage + binary_damage_ratio")
    )	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[2], "~ . - ratio_damage + binary_damage_ratio"))	%>%
    gam(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[4], "~ . - ratio_damage + binary_damage_ratio"))	%>%
    gam(data = combined_tbl, family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE,
  keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------


## ----------------------------------------------------------------------
##
## Table A5
## Result IV: >70 residential ratio subsample
##
## ----------------------------------------------------------------------
## Panel A: Damage ratio
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[3], " ~ ."))	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ .")
    )	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[2], " ~ ."))	%>%
    # update(str_c(depvar_vec_combined[1], " ~ ."))	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7)),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[4], " ~ ."))	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table A5: Panel A
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE,
  keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------
## Panel B: Dummy regressions
## ----------------------------------------------------------------------
rslt_lst_dummy = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[2], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7)),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[4], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(ratio_residential >= 0.7), family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table A5: Panel B
stargazer(
  rslt_lst_dummy,
  type = "text",
  align = TRUE,
  keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------


## ----------------------------------------------------------------------
##
## Table A6
## Result III: Close-palace subsample (PalaceDist <= 10 km)
##
## ----------------------------------------------------------------------
## Panel A: Damage ratio
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl	%>%	filter(palace_dist <= 10), keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[3], " ~ ."))	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ .")
    )	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[2], " ~ ."))	%>%
    # update(str_c(depvar_vec_combined[1], " ~ ."))	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10)),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[4], " ~ ."))	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10), family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table A6: Panel A
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)

## ----------------------------------------------------------------------
## Panel B: Dummy regressions
## ----------------------------------------------------------------------
rslt_lst_dummy = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl	%>%	filter(palace_dist <= 10), keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10), family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[2], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10)),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[4], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3])
    )	%>%
    gam(data = combined_tbl	%>%	filter(palace_dist <= 10), family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table A6: Panel B
stargazer(
  rslt_lst_dummy,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------

## ----------------------------------------------------------------------
##
## Table A7
## Result VI: with spatially lagged Damage
##
## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . + LaggedDamage")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[3], "~ . + LaggedDamage"))	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ . + LaggedDamage")
    )	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[2], "~ . + LaggedDamage"))	%>%
    gam(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[4], "~ . + LaggedDamage"))	%>%
    gam(data = combined_tbl, family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table A7: Panel A
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE,
  keep = c(treatment_vec, robust_treatment_vec, "LaggedDamage")
)

## ----------------------------------------------------------------------
## Panel B: Dummy regressions
## ----------------------------------------------------------------------
rslt_lst_dummy = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3], " + LaggedDamage")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3], " + LaggedDamage")
    )	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[3], "_1993 ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3], " + LaggedDamage")
    )	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[2], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3], " + LaggedDamage")
    )	%>%
    gam(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%
    update(
      str_c(depvar_vec_combined[4], " ~ . - ", treatment_vec[1], " + ", robust_treatment_vec[3], " + LaggedDamage")
    )	%>%
    gam(data = combined_tbl, family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table A7: Panel B
stargazer(
  rslt_lst_dummy,
  type = "text",
  align = TRUE,
  keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------


## ----------------------------------------------------------------------
##
## Table A9
## Result I: ANA result---continuous damage, 5th poly
##
## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[4]]	%>%	update(str_c(depvar_vec_combined[3], " ~ ."))	%>%
    glm(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[4]]	%>%	update(str_c(depvar_vec_combined[3], "_1993 ~ ."))	%>%
    glm(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[4]]	%>%	update(str_c(depvar_vec_combined[2], " ~ ."))	%>%
    # update(str_c(depvar_vec_combined[1], " ~ ."))	%>%
    glm(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[4]]	%>%	update(str_c(depvar_vec_combined[4], " ~ ."))	%>%
    glm(data = combined_tbl, family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
## Table 1: Panel A
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)

## ----------------------------------------------------------------------
##
## Table A10
## Result XI: Baseline regressions without some of AP and IP polynomials
##
## ----------------------------------------------------------------------

## --- Polynomials and splines only IP
distance_term_specification_IP = list(
  str_c("poly(", distance_termz_base[1], ", degree = ", polynomial_degree, ", simple = TRUE, raw = TRUE)"),
  str_c("s(", distance_termz_base[1], ")")
)

## --- Polynomials and splines only AP
distance_term_specification_AP = list(
  str_c("poly(", distance_termz_base[2], ", degree = ", polynomial_degree, ", simple = TRUE, raw = TRUE)"),
  str_c("s(", distance_termz_base[2], ")")
)


mdl_lst = c(mdl_lst,
            ## ------------------------------------------
            ## Additional models for this table
            ## ------------------------------------------
            ## Without polynomials of longitude & latitude
            polynomial_base_wo_lonlat = chr2fml_felm(
              "outcome",
              idv_list = list(
                treatment_vec[1],
                spatial_covariates,
                distance_term_specification[[1]]
              ),
              fe_list = fe_term,
              se_cluster = "row_count"
            ),
            ## Without polynomials of Aiming Points
            polynomial_base_wo_AP = chr2fml_felm(
              "outcome",
              idv_list = list(
                treatment_vec[1],
                spatial_covariates,
                distance_term_specification_IP[[1]],
                spatial_term_specification[[1]]
              ),
              fe_list = fe_term,
              se_cluster = "row_count"
            ),
            ## Without polynomials of Imperial Palace
            polynomial_base_wo_IP = chr2fml_felm(
              "outcome",
              idv_list = list(
                treatment_vec[1],
                spatial_covariates,
                distance_term_specification_AP[[1]],
                spatial_term_specification[[1]]
              ),
              fe_list = fe_term,
              se_cluster = "row_count"
            ),
            ## Without polynomials of AP and IP
            polynomial_base_wo_APIP = chr2fml_felm(
              "outcome",
              idv_list = list(
                treatment_vec[1],
                spatial_covariates,
                spatial_term_specification[[1]]
              ),
              fe_list = fe_term,
              se_cluster = "row_count"
            ),
            ## Without polynomials of lon&lat, AP, IP
            polynomial_base_wo_LonLatAPIP = chr2fml_felm(
              "outcome",
              idv_list = list(
                treatment_vec[1],
                spatial_covariates
              ),
              fe_list = fe_term,
              se_cluster = "row_count"
            )
)

## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## 1b. DV = dummy_ninka, LPM
  ## wo_LonLat
  Dummy_Ninka_LPM = mdl_lst[[5]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## 1c. DV = dummy_ninka, LPM
  ## wo_AP
  Dummy_Ninka_LPM = mdl_lst[[6]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## 1d. DV = dummy_ninka, LPM
  ## wo_IP
  Dummy_Ninka_LPM = mdl_lst[[7]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## 1e. DV = dummy_ninka, LPM
  ## wo_APIP
  Dummy_Ninka_LPM = mdl_lst[[8]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## 1f. DV = dummy_ninka, LPM
  ## wo_LonLatAPIP
  Dummy_Ninka_LPM = mdl_lst[[9]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE)
)

## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE,
  keep = c(treatment_vec, robust_treatment_vec)
)

## ----------------------------------------------------------------------
##
## Table A.11
## Result II: Regression with 5th ordered polynomials of prewar covariates
##
## ----------------------------------------------------------------------
## Result list object
## ----------------------------------------------------------------------
rslt_lst = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ . + ", extra_covariates_specification[2])	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = combined_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[3], " ~ . + ", extra_covariates_specification[2]))	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[3], "_1993 ~ . + ", extra_covariates_specification[2]))	%>%
    gam(data = combined_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst[[3]]	%>%
    update(str_c(depvar_vec_combined[2], " ~ . + ", extra_covariates_specification[2]))	%>%
    # update(str_c(depvar_vec_combined[1], " ~ . + ", extra_covariates_specification[2]))	%>%
    gam(data = combined_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst[[3]]	%>%	update(str_c(depvar_vec_combined[4], " ~ . + ", extra_covariates_specification[2]))	%>%
    gam(data = combined_tbl, family = count_link)
)

## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
stargazer(
  rslt_lst,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)
## ----------------------------------------------------------------------


## ----------------------------------------------------------------------
##
## Table A12
## Result VIII: Baseline regressions with zip code-level data
##
## ----------------------------------------------------------------------
## Baseline formula objects
## ----------------------------------------------------------------------
mdl_lst_zipcode = mdl_lst

## ----------------------------------------------------------------------
## Panel A: Damage ratio
## ----------------------------------------------------------------------
rslt_lst_zipcode = list(
  ## ----------------------------------------------------------------------
  ## 1a. DV = dummy_ninka, LPM
  Dummy_Ninka_LPM = mdl_lst_zipcode[[2]]	%>%
    update(
      str_c(depvar_vec_combined[3], " ~ .")	%>%	as.formula()	## formula obj to feed felm()
    )	%>%
    felm(data = zipcode_tbl, keepCX = TRUE, psdef = FALSE),
  ## ----------------------------------------------------------------------
  ## 1b: DV = dummy_ninka, probit/logit link
  Dummy_Ninka_GAM = mdl_lst_zipcode[[3]]	%>%	update(str_c(depvar_vec_combined[3], " ~ ."))	%>%
    gam(data = zipcode_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 1c: DV = dummy_ninka (post-1993 period), probit/logit link
  Dummy_Ninka_GAM_post1993 = mdl_lst_zipcode[[3]]	%>%	update(str_c(depvar_vec_combined[3], "_1993 ~ ."))	%>%
    gam(data = zipcode_tbl, family = binom_link),
  ## ----------------------------------------------------------------------
  ## 2. DV = log(ninka_days), Gaussian link
  NinkaDays_Linear = mdl_lst_zipcode[[3]]	%>%	update(str_c(depvar_vec_combined[2], " ~ ."))	%>%
    # update(str_c(depvar_vec_combined[1], " ~ ."))	%>%
    gam(data = zipcode_tbl),
  ## ----------------------------------------------------------------------
  ## 3. DV = n_ninka, (quasi-)Poisson link
  N_Ninka_GAM = mdl_lst_zipcode[[3]]	%>%	update(str_c(depvar_vec_combined[4], " ~ ."))	%>%
    gam(data = zipcode_tbl, family = count_link)
  ## ----------------------------------------------------------------------
)
## ----------------------------------------------------------------------
## Print regression estimates
## ----------------------------------------------------------------------
stargazer(
  rslt_lst_zipcode,
  type = "text",
  align = TRUE, keep = c(treatment_vec, robust_treatment_vec)
)




## EOF



